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5> ■ ABSTRACT 

d ; 

We present a detailed study of the critical properties of the 2-D XY model with 
maximal frustration in a square lattice. We use extensive Monte Carlo simulations 
q \ to study the thermodynamics of the spin and chiral degrees of freedom, concen- 
• rH ! trating on their correlation functions. The gauge invariant spin-spin correlation 

X: 

functions are calculated close to the critical point for lattice sizes up to 240 x 240; 
the chiral correlation functions are studied on lattices up to 96 x 96. We find that 
the critical exponents of the spin phase transition are v = 0.3069, and r\ = 0.1915, 
which are to be compared with the unfrustrated XY model exponents v = 1/2 
and 7] = 0.25. We also find that the critical exponents of the chiral transition 
are v x = 0.875, 2(3 = 0.1936, 27 = 1.82, and 27 / = 1.025, which are different 
from the expected 2-D Ising critical exponents. The spin-phase transition occurs 
at Tf/m = 0.446 which is about 7% above the estimated chiral critical temperature 
Tz 2 = 0.4206. However, because of the size of the statistical errors, it is difficult 



to decide with certainty whether the transitions occur at the same or at slightly 
different temperatures. Finally, the jump in the helicity modulus in the fully frus- 
trated system is found to be about 23% below the unfrustrated universal value. 
The most important consequence of these results is that the fully frustrated XY 
model appears to be in a novel universality class. Recent successful comparisons 
of some of these results with experimental data are also briefly discussed. 

PACS numbers: 5.70.Jk, 74.50. +r, 64.60Fr. 



2 



1. Introduction. 



The critical behavior of the uniformly frustrated 2-D XY model has been stud- 
ied extensively in recent years, both theoretically 1-25 and experimentally 26-32 . 
This theoretical interest has been due to the rich variety of possible novel critical 
phenomena that can appear in this model depending on the frustration parameter 
/ = v/Vi with p and q relative primes. Experimentally an understanding of the 
phase transition(s) that occurs in this model is important to describe the physics of 
two-dimensional periodic arrays of Josephson junctions, 26-32 and two-dimensional 
superconducting wire networks 33 " 34 , both in the presence of frustration / = $/$o- 
Here <E> is the average flux per plaquette normalized to the superconducting quan- 
tum of flux $o — h/2e. These arrays can be manufactured with high precision using 
modern photolithographic techniques. Of particular interest is the / = 1/2 fully 
frustrated 2-D XY model (FFXYM). This model has a continuous U(l) abelian 
symmetry, and a discrete Zi symmetry leading to the possibility of true long-range 
order in two dimensions. In contrast, the unfrustrated 2-D XY model (XYM) only 
possesses a continuous U(l) abelian symmetry: its low temperature phase is char- 
acterized by quasi- long range order rather than true long-range order 35-38 . In spite 
of the many experimental and theoretical studies of the FFXYM, there are several 
questions that remain to be resolved. For example, it is not clear whether one 
phase transition exists at the critical temperature T c which is a combination of a 
Berezinskii-Kosterlitz-Thouless (BKT) type transition for the U(l) symmetry plus 
an Ising-like transition for the Zi symmetry, or whether there are two successive 
phase transitions at critical temperatures an d Tz 2 - Even the order in which 

they may occur is controversial. More importantly the nature of the transitions, 
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as characterized by their critical properties, is not yet fully understood. 

In their original work Teitel and Jayaprakash suggested 4 that in a square lattice 
the two transitions occurred very close in temperature. They carried out Monte 
Carlo (MC) simulations to calculate the helicity modulus T and the specific heat 
C as a function of temperature and lattice sizes L x L, with L up to L = 32. 
They found that the maximum of the specific heat appeared to increase as lnL, 
characteristic of a 2-D Ising-like transition. Related studies in the triangular lattice 
antiferromagnetic XYM, which could be expected to be in the same universality 
class as the FFXYM, indicate that there is a combination of BKT and Ising-like 
transitions 6 . In ref. 6, the two transitions appear to take place at the same T c 
while in ref. 7 they are within 2% of each other, with Tz 2 > Tjjny In another 
investigation Berge et al. 14: introduced a frustrated XYM with variable frustration 
on a square lattice. In this model the couplings along the columns are chosen with 
strength J, while those along every other row have strength —fiJ, with < /i < 1. 
From a MC analysis of the specific heat they surmised that for fi < 1 the model 
has separate Ising and BKT ordering with Tz 2 < while for the FFXYM 

(/i = 1) the two transitions appear to merge into one. The Berge et al. model was 
studied in detail by Eikmans et al. 16 who carried out MC calculations of the helicity 
modulus, which is more sensitive to possible BKT-like ordering, and interpreted 
their results using a Coulomb gas picture. In another thermodynamic MC study 
of the related uniformly frustrated square lattice 2-D Coulomb gas model, Grest 17 
carried out simulations for frustrations / = 0, 1/2, 1/3 and 1/4 and with L up to 
50. He found that in the fully frustrated case the jump in the inverse dielectric 
constant e" 1 is different from the XYM case. Specifically, the jump in e^ 1 occurs 
at T CG = 0.129 ± 0.002 and takes the value e" 1 = 0.63 ± 0.03, which is larger than 
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the XYM universal value of 0.52, and agrees with Minnhagen's conjecture. The 
determination of the jump was based, however, on the criterion used for the XYM. 
Grest found, as in previous studies, that the specific heat grows logarithmically 
with L. It is significant that his results appear to indicate a clear separation of the 
two critical temperatures with Tz 2 > ^t/(i), contrary to previous conjectures. 

Granato 19 et al. have studied the Z2 critical behavior of a coupled XY-Ising 
system using MC and MC-transfer matrix calculations. An important finding in 
this study is the chiral critical exponent v x ~ 0.85(3), which is clearly different 
from the 2-D Ising model value of v = l. 18 Furthermore, they found that the XY 
and Ising transitions occur at essentially the same temperature. Lee et al. 20 carried 
out MC simulations of the FFXYM in the square and triangular lattices and found 
that v is also different from the 2-D Ising model result. 

In the XYM the nature of the BKT phase is characterized by the approximate 
analytic expression for the spin-spin correlation functions 36 ' 37 . However, unlike in 
the XYM case, it has proven to be very difficult to calculate the correlation func- 
tions for the FFYXM analytically. This difficulty exists partly because in order to 
carry out the calculations one needs to include the basic excitations of the frus- 
trated problem, which consist of different types of fractional charges as well as the 
Ising model related domain walls 11 . Nonetheless, it has been possible to extract 
some qualitative information about the critical properties using techniques such as 
the renormalization group approximation 9 applied to an effective hamiltonian ob- 
tained from a Hubbard- Stratonovich transformation of the FFXYM 8 or by general 
symmetry arguments 10 . One worry about the effective hamiltonian is that it does 
not explicitly contain the same elementary excitations as the original FFXYM, 
such as the fractional charges. All of the studies mentioned above have mostly 
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concentrated on calculating thermodynamic quantities, for it has been difficult to 
separate the Z2 from the U(l) contributions. 

The purpose of this paper is to fill this gap by explicitly calculating the U(l) 
and Z2 correlation functions as well as the separate Z2 contribution to the magnetic 
properties. We should mention at the outset that these calculations are significantly 
more demanding than the thermodynamic calculations and are now possible be- 
cause of improved algorithms and computer power. One further complication is 
that at present there is no available analytic theory for the / 7^ case that could 
suggest what form these correlation functions should have and we need to make an 
ansatz for them. Generally, we can either assume that they decay exponentially 
or algebraically with distance. We use different statistical measures to test for 
the two possibilities. If our MC results for the correlation functions are consistent 
with an exponential decay we extract a correlation length £(T), while if they are 
consistent with a power law decay we extract the corresponding i](T) exponent. 
In the case that £(T) diverges at T c from above it can diverge as a power law or 
with the BKT form ~ exp {B(T — T c ) -Iy ). In the / = case the critical exponent 
v(f = 0) = 1/2. 36 ' 37 In the low temperature phase of the XYM the correlation 
function decays algebraically with distance r as ~ r~' n , where the rj exponent is a 
continuous function of T and takes the universal value rj(f = 0, Tbkt) = 1/4. Sev- 
eral experiments have confirmed the / = picture and the values of the measured 
critical exponents 40-42 agree well with those predicted by theory. In addition, 
recent MC simulations have provided an accurate evaluation of the / = XYM 
critical exponents 43-47 . The most recent 46 high statistics estimates for / = are 
v = 0.4695(1) and i] = 0.235, with the critical temperature Tbkt — 0.8953. 

In order to understand the nature of the phase transitions in the FFXYM 
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we have studied a variety of quantities, several of which separately describe each 
particular symmetry. The thermodynamic quantities calculated are: the helicity 
modulus, T, and the square of both the staggered chiral magnetization, A4 2 ., and 
susceptibility, x\- We have carried out an extensive analysis of the gauge invariant 
U(l) correlation function, 9(u(i) ( r ) an d their corresponding even and odd coherence 
lengths (to be defined below). These calculations have allowed us to extract the 
U(l) critical temperature, 2|m)> and its critical exponents v and 77. For the Z2 
freedoms we calculated the chiral correlation function, g x (r), and its corresponding 
coherence length, £ x , which allowed us to estimate the critical exponent v x and the 
critical temperature Tz 2 ■ Our result for the exponent v x is in very good agreement 
with the recent MC transfer matrix calculation 25 . 

We will now outline the main results of our study. Our extensive analysis is 
consistent with a U(l) BKT-type transition but with exponents z/(f=l/2)= 0.3069 
and 77 (f=l/2, T c ) = 0.1915. These results clearly differ from those obtained in 
the XYM case 36 > 37 > 43 . We have also calculated the Z2 critical exponent 2(3 = 
0.1936(35) for M 2 S , 2 7 = 1.82(13) and 2 7 / = 1.025(79) for x % and the coherence 
length exponent v x = 0.875. These exponents are also different from those expected 
for a 2-D Ising model. The critical temperatures found in our study are 71/(1) = 
0.446 and Tz 2 = 0.4206. One could be tempted to say that the transitions take 
place at two different temperatures, and this may indeed be the case. However, 
after a detailed assessment of the size of the statistical errors from the nonlinear fits 
and considering the small difference between the two temperatures we can not be 
certain if they are different or not. Furthermore, the transitions are reversed from 
their expected order. We suspect that more extensive simulations with algorithm 
improvements, better statistics and larger system sizes are needed to clarify this 
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point. The results mentioned above were obtained from extensive MC simulations 
on a square lattice of size L, with L ranging from L = 8 up to L = 240, and 
with periodic boundary conditions. What emerges from our results is that the 
FFXYM is in a novel universality class different from either a pure XY or an Ising 
universality class. A brief description of some of the results presented here has 
appeared elsewhere 22 . 

We have also recently reanalyzed the experimental results of the Delft group 48 
for / = and / = 1/2. We have concluded that the values of rj(f = 0) = 1/4 and 
i](f = 1/2) = 0.1915 ~ 1/5 are in good agreement with the experimental data. 
However, the fits of the experimental resistance versus temperature data can not 
distinguish between a v(f = 1/2) = 1/3 from a v(f = 1/2) = 1/2. Moreover, as 
mentioned above, recent MC-transfer matrix work has provided further evidence 
that the chiral exponents are not equal to the 2-D Ising model exponents and even 
quantitatively the value of v x has begun to converge on values close to 0.85. 

The organization of this paper is as follows: In section 2.1 we define and 
briefly review the general properties of the uniformly frustrated 2-D XYM. In 
section 2.2 we define the thermodynamic quantities calculated in this paper while 
in section 2.3 we give the expressions for the calculated gauge invariant U(l) and Z2 
zero-momentum correlation functions, of central interest here, together with their 
possible asymptotic behaviors. In Section 3.1 we describe briefly the MC algorithm 
used in our calculations. Since there are no analytic results for the correlation 
functions to guide our analysis, we proceed by developing an approach that consists 
of using several independent checks of the results obtained. As a test, in subsection 
3.2, we successfully apply our strategy to the unfrustrated XYM and compare our 
results to those obtained in the more extensive recent MC studies 43-46 . In section 
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4 we present the bulk of our numerical MC results applied to the FFXYM. In 4.1 
we discuss the thermodynamic results for both the U (1) and Z2 freedoms. In 4.2(a) 
we give the correlation function results for the U(l) freedoms, including a finite 
size scaling analysis for the correlation length. In section 4.3(b) we present the 
corresponding correlation function results for the Z2 freedoms. Finally in section 

5 we present a critique of our results and a possible outlook for the future. 

2. The fully frustrated XY model 

2.1. Definition of the model. 

The uniformly frustrated 2-D XYM is defined by the hamiltonian 

H = - Jcos(e(r)-e(f/)+f(r,f/)y (1) 

<f,f/> 

where 6(f) is the angle at site f, < f,ff > stands for a sum over nearest- neighbor 
lattice sites, and J is the exchange constant. In the Josephson junction array 
representation of the model in a transverse magnetic field, the bond variables 
/(r, ft) are given by the line integral f(f, f /) = |^ JZ 1 A - dl, with A the magnetic 
vector potential. For uniform frustration these bond angles are required to satisfy 

E f^ = T j A-dl = 2Kf. (2) 

Plaquette ° plaquette 

The hamiltonian defined in Eq. (1) is invariant under the transformation, 9(f) — > 
9(f)+2nn(f) and f(f, ff) — > f(f, ff)+2n[n(rf)— n(f)}, where n(f) and n(ff) are inte- 
ger numbers. 
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Choosing the gauge A = (—By, 0, 0) so that B = Bz and assuming a square 



lattice, the bond angles f(r,rt) are given by 



/(f, ft) = T^fU + for f = r ± a i and 
f(f,ff) = for r = f ± aoj. 



(3a) 



(36) 



Here / = Ba^/^ with r = (ia ,ja ), i,j integers and a the lattice spacing. The 
uniformly frustrated model is periodic in / with period one, and with reflection 
symmetry about / = 1/2. The XYM corresponds to the unfrustrated / = case. 

The fully frustrated case corresponds to / = 1/2. The effect of / in this case 
is to produce alternate rows with ferro- and antiferromagnetic couplings, while the 
couplings along the columns are all ferromagnetic. Each plaquette has one antifer- 
romagnetic and three ferromagnetic bonds, or vice versa, leading to a ground state 
that has a two-fold degeneracy with half-integer vortices of opposite circulation or 
chirality. 1 Thus, the system displays two symmetries: the underlying continuous 
U(l) abelian symmetry for the phases and a discrete Z2 or Ising-like symmetry 
associated with the chiral degrees of freedom. 

2.2. Thermodynamic properties 

The helicity modulus T is defined by the response of the system to a twist in 
the spins at its boundaries. In our case T is calculated explicitly from the formula, 
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+ k^[\ £ ^sin(^(r)-^f/)+/(r,f/) 

\<r,r/> 

where ( ) stands for a thermal average, and ks is Boltzmann's constant, and x? p = 
Another informative quantity associated with the supercurrent loops 
around a plaquette in a Josephson array is the staggered magnetization 



M. 



stagg. 



N 



Rx~\~Ru 



^V(R) 



J2 mn(e(f)-9{f/)+f(f,rt)) ), (5) 

<r,f f>€V{R) 



where (it^, give the coordinates at the center of the plaquette V . The index 

— * 

V(R) runs from 1 up to N, the total number of plaquettes in the lattice. 

We now turn to the definition of the quantities associated with the chiral 
degrees of freedom. The chirality of a plaquette gives the direction of circulation 
of the supercurrents induced by frustration. Each plaquette has a definite chirality 
which can be ±1, and it is calculated from 1 



x(R) = sign sin^(r) — 6(f t) + f{r,ff)j 

<r,f/>eV(R) 



(6) 



with the dual lattice vector R = [(i + l/2)a , (j + l/2)a ] with i,j integer numbers. 
At zero temperature the chiralities are ordered like a 2-D Ising antiferromagnet. 
At finite temperatures there are line or domain wall defects separating regions with 
different chiralities. 

The order parameter describing the Z2 phase transition is the staggered chiral 
magnetization, defined by 



(7) 



It is difficult to study this quantity numerically since it oscillates rapidly between 
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positive and negative values. A more stable quantity to study is Binder's second 
order cumulant, M. 2 S and its fluctuations. 7 These fluctuations are given by the 
square of the staggered chiral susceptibility , calculated as 



Our numerical results for the thermodynamic quantities characterizing the chi- 
ral degrees of freedom will be presented in section 4.1(b). It will be seen that close 
to the critical region their behavior is quantitatively different from an Ising ferro- 
magnet on a square lattice. 

2.3. Correlation functions 

It is known that the hallmarks of the BKT ordering can be given in terms of 
the phase correlation functions. An analytic evaluation of these quantities appears 
to be mathematically intractable for uniform frustration. Nonetheless, for random 
frustration it has been possible to analytically calculate the correlation functions 
at low temperatures in the limits where the density x / of frustrated plaquettes 3 is 

< 1 or ~ 5. Given that the hamiltonian of the frustrated XYM is gauge 
invariant, the phase correlation functions should also be gauge invariant. The 
correlation functions are defined along a path connecting the correlated spins and 
are therefore path-dependent. In fact, gauge invariant correlation functions along 
two different paths differ by the total amount of frustration enclosed by the two 
paths. The gauge invariant phase correlation function along a path T that accounts 
for the frustration in the system is given by 2 ' 3 




(8) 
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In Fig. 1 we show two possible trajectories r^(r, ft) and r#(r*, ft) joining the points 
f and ft. The phase introduced in Eq.(9) when going around the trajectory A is 
Y[r A e*-^ s ' s/ ) = exp^E r ^ f(s, st)j, while if one follows trajectory B the phase fac- 
tor is given by 
Ilr s e*-^ s ' s/ ) = exp^E rs f(s, s The total frustration contained in the area 
enclosed by both trajectories is then 

f(s,st)+ f(s,st) = 27rMf, (10) 

<s,s/>er A <s,s/>er B 

where M represents the number of elementary plaquettes inside the area encircled 
by the paths, and / is the frustration of each plaquette. In this case the phase 
in Eq.(9) is shifted by an additional amount 2tt Mf when the correlation along 
trajectory B is calculated instead of A. 

To evaluate the effect of frustration on the correlation functions at low tem- 
peratures one can use duality transformations to obtain a lattice Coulomb gas 
representation of the model. After doing so one obtains 3 

9u(i)(r,f) =esxp^i^2^2-n(f)e(r-R)f(R) x g c . g Xf,ff), (11) 

r R 

where g c . g .(^ ?') l $ the lattice Coulomb gas correlation function for the XYM 37 and 

— * — * 

f(R) is the frustration at the plaquette with center at the dual lattice site R. The 
number n(f) is zero everywhere, except at f and ft, where it takes the values n{f) = 
—n(ft) = l. 3 The angular potential 6(.R) is given by 37 iO(z) = \n(z) - G(\z\), 
for large R. Here z = R x + iR y and R = (R x ,R y ). Notice that in Eq.(ll) 
there is an extra phase factor appearing in the phase correlation function. This 
factor weights the contributions to the correlation function coming from different 
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trajectories going between r and r', and it is a consequence of gauge invariance or, 
equivalently, the Aharanov-Bohm effect. We shall see that this extra phase factor 
in the correlation function appears naturally in the results discussed in Section 
4.1(a). 

In the evaluation of the correlation functions we can have important contribu- 
tions from more than one Lyapunov exponent, which makes the extraction of the 
largest exponent difficult. However, this problem is not present if we evaluate the 
zero momentum correlation function defined by 43 

9o (r) =< S av (i) ■ S av {i + r) >, (12) 

where 

§av(i) = ^Y,§(i,j) (13) 

Ly 3 

is the average spin along the ith-column. In this case r denotes the distance 
between the columns being correlated. 

From the definition of the gauge- invariant phase correlation function, Eq. (9), 
the expression for the zero momentum correlation function is 

/ 1 Ly Lx 1 \ 

9u(l)(r) = ( j-j- cos ( i+fJ " Oij + *ti + ~) r ) ) , ( 14 ) 

where we have used Eqs.(3). 

A similar reasoning applies to the zero momentum chiral correlation function 
given by 

/ 1 L x -lL y -l \ 

9x ^ = \ TT ^ ^ Xi+rJ XhJ ) ' ^ 15 ^ 
\ x y i=i j=i I 

We have used equations (14) and (15) to evaluate the correlation functions in our 
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numerical simulation. 

As mentioned before there are no known explicit analytic expressions for <7t/(i)(r)| 
and g x ( r )- Nonetheless, we can make an ansatz for the analytic form of these corre- 
lations close to the critical region. Our ansatz is based on what is known about the 
XYM and about general properties of standard second order phase transitions. In 
the XYM as T — > T^ KT the asymptotic form of the spin-spin correlation function 
for r ^> 1 is 

with r/Q = 1/4 and the coherence length diverging as 

In the XYM the critical exponent vq = 1/2. 36 In the low temperature phase (T < 
Tbkt) the long distance correlation function decays as 

»M = ;§r < 18 > 

The exponent 770 is a function of temperature, representing a continuous line of 
critical points. In the XYM 770 takes the universal value t)q{T = Tbkt) = 1/4. 36 ' 37 
This result is directly related to the universal jump predicted for the superfluid 
density. 38 

In the disordered phase of the 2-D Ising model (T > Tj), the asymptotic 
behavior of the correlation function close to Tj is given by 49 

9i(r) = ^ x (19) 
15 



again with 77/ = 1/4 and with a power law divergence for the coherence length 

&C0 = (jr^rp (20) 

where the critical exponent z/j = 1. In the ordered phase (T < Tj), the asymptotic 
behavior of the correlation functions is given by 49 

9i(r) = z~ fl + < M i > 2 , (21) 

valid for ej = Tl j7 r << 1 and ej r < 1. The correlation function critical exponent 
at Tj is 771 = 1/4, as in the XYM. 

Based on previous studies of the thermodynamics of the FFXYM it is reason- 
able to assume that their asymptotic behavior for / = 1/2 can be described by 
either BKT or Ising-like forms described above. In our calculations we checked for 
the best fits to our MC data by either form. 

Having discussed the expected analytic forms for the different correlation func- 
tions of interest, let us now turn to a discussion of their numerical evaluation. We 
should notice first that, strictly speaking, for finite lattices the asymptotic behav- 
ior of the correlation functions is not accessible. Even in rather large lattices the 
subleading power law behavior of the correlation functions can be non-negligible. 
Thus, the evaluation of the coherence length extracted from a numerical calculation 
of the correlation function is nontrivial. A common procedure is to take periodic 
boundary conditions and then fit the behavior of the correlations to 

G(r)=g(r)+g(L-r), (22) 

where g(r) is any of the correlation functions of interest. We should note that in 
general it is not sufficient to account for the closest images to the source along the 
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r direction but we must also account for the images farther away as well as those in 
the direction transverse to r. In fact, their contribution becomes more important 
as the coherence length £ grows since the number of relevant images increases. 

3. Calculational strategy and test 

Our strategy is to carry out several independent consistency checks of our 
results, for there are no analytic results with which to guide the analysis. To test the 
reliability of our consistency checks, we start by applying them to the extensively 
studied XYM. Although there is a basic consensus about the physical nature of 
the BKT transition, relatively reliable and thorough nonperturbative numerical 
studies of the critical exponents of the XYM became available just recently. 43-47 
Here we tried to follow some of the basic ideas of these approaches, in particular 
the one used above T c , complemented with other tests implemented here. We must 
stress that we are on less firm ground in the FFXYM case than in the XYM and 
thus we need extra consistency checks that were not needed in the XYM studies. 

3.1. MC algorithm 

Different acceleration algorithms that have worked out well in the XYM were 
constructed with special regard to the nature of the basic excitations in the model. 
In the FFXYM we have a less definitive idea about the basic excitations in the 
model and therefore the same type of algorithms have not proven any more efficient 
than the standard Metropolis approach. 49 ' 50 Our simulations were then carried 
out using the standard Metropolis algorithm in square lattices of sizes L x L with 
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L = 8, 16, 32, 60, 72, 84, 96, 180, and 240, with periodic boundary conditions. The 
lattice sizes L = 120, and L = 180, 240 were considered only at two temperatures 
that are about 3% and 2% away from the estimated critical point, respectively. 
These relatively large lattices were studied in order to get a better estimate of the 
U(l) correlation function critical exponents. For the thermodynamic and chiral 
degrees of freedom a full set of temperature values was considered for L up to 
L = 84. Since we expect to have critical slowing down similar to that seen in the 
XYM, which has decorrelation times growing as £ 2 , we performed reasonably long 
runs, although no detailed attempt was made to calculate the FFXYM dynamic 
critical exponent. Nonetheless, our consistency and self-consistency checks give 
support to the reliability of our results. The equilibration time of a typical run 
was at least 10K MCS/angle far from criticality and at least twice as much for 
temperatures close to T c . The statistics were calculated from runs of at least 50K 
MCS/angle, and up to 290K MCS/angle. Details of the length of the runs are 
given in the tables. It is important to note here that since we are interested in 
extracting critical exponents from nonlinear fits, plots of the results are sometimes 
not as informative as looking at the numbers themselves. 



3.2. Unfrustrated 2-D XY model 



We begin by presenting our results for the XYM together with the tests of 
our consistency checks for both T° and the correlation functions. We compare 
our XYM results to those obtained recently by more extensive analysis 46 . In the 
next section we shall present the bulk of the results of our calculations with the 
FFXYM. Note that in terms of the correlation functions the difficult calculations 
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are those for the U(l) symmetry, for in that case the correlation length may diverge 
exponentially rather than algebraically as one gets close to T c . 

We now describe the numerical approach to calculate the coherence length, the 
critical temperature and critical exponents from the U(l) correlation functions. 
We have to some extent repeated the recent, more extensive calculations for the 
XYM 43-46 critical exponents for T > Tbkt, and we have extended the calculations 
to the T < Tbkt region. To facilitate the comparison of our results, given in 
Table 1A, to previous findings we have summarized the results of references 43- 
46 in Tables IB and 1C. Basically, we have followed the method of analysis used 
in those references, although the lattices we have simulated are not as large as 
the ones considered there. However, to reduce the finite size effects and to insure 
meaningful results for the correlation functions, we kept the ratio L/^q > 4 in all the 
calculations. We note that, even though the critical temperatures in the FFXYM 
is about Tbkt/2, we kept this ratio at L/£ > 5. The calculational procedure is the 
following: First, the periodic form of the zero- momentum correlation function go(r) 
given in Eqs. (14) and (22), with / = 0, was calculated in the high temperature 
phase. Next, we carried out unconstrained 3-parameter nonlinear fits of the data to 
the form given in Eq(16). From these fits we determined £,o(T) and the parameters 
i]q and Aq. To further check the consistency of the nonlinear fits we performed 
linear fits to the MC data of the form 



varying the values of t]q until we reached a minium for the x 2 function. We found 
that £o(T) is systematically above the values of Ref. 43 (for the comparison see 
Fig. 2), and that rjo(T) oscillates nonmonotonically close to the critical point, mak- 
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(23) 



ing its determination difficult. We also carried out an unconstrained 4-parameter 
nonlinear fit to the data to extract Tbkt and vq, assuming a BKT form for £o(T). 
As pointed out in a recent finite size scaling analysis 45 of the XYM susceptibility, 
it is very difficult to distinguish between a BKT form and a power law divergence 
on the basis of data obtained from MC simulations on lattices up to L=256, even 
for temperatures 1% away from the critical point. Thus, we carried out an ad- 
ditional unconstrained 3-parameter nonlinear fit to a power law form [Eq. (20)]. 
The results of these fits are summarized in the first and second columns of Table 
1A. Before proceeding with a discussion of these results, let us emphasize that we 
carried out a careful analysis of the stability of the parameters obtained from the 
fits by trying to make sure that the values obtained correspond to the minimum 
of the x 2 function, within the statistical errors of our simulations. More details of 
the fitting procedure and related analyses will be discussed below. 

Let us now turn to the results obtained by assuming a BKT form for £o(T). 
We find that the values of the parameters Aq and -Bo in the first column of Tables 
1A and IB agree well, while the values of vq and Tbkt differ by 4% and 0.75%, 
respectively. The results obtained assuming a power law form are given in the 
second column of Tables 1A and IB. The values of T c and z/ Q are within 1% and 
9%, respectively, whereas the values for the A^s are the same. All of these results 
indicate that there is good agreement between our results and those obtained from 
the more extensive MC simulations, in spite of the fact that our simulations are 
in smaller systems and have less statistics. It is important to realize that fitting 
the coherence length to a BKT or to a power law form leads to x 2 /dof- values of 
comparable quality. Thus, from this analysis alone one cannot decide which of 
the two fits is the correct one. To sort out this problem we also have calculated 
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the low temperature go(r,T). The correlations, for L = 60 lattices, were fitted to 
go(r) ~ r~ Vo and 



9o(r) 



Vo p a r I i t r \-Vo p ao{L-r) 



(24) 



at 9 different temperatures. The results of these fits are given in Tables 2A and 
2C. It is found that the exponential fits appear to be of better quality, and the 
corresponding values of rjo(T) are systematically above those obtained from the 
algebraic fits. However, the important point is that the values of the exponents 



trusted. Moreover, the algebraic fits agree with the low T spin wave prediction, 
r)o(T) ~ T. In Fig. 3 (a) we show the results for 770 as a function of temperature. To 
further analyze the nature of the low temperature phase, we calculated t]q(Tbkt) = 
r/oc as a function of lattice size, using the Tbkt obtained at high temperatures and 
for reasonably long runs. On the other hand, the exponential fits yielded values of 
r/oc systematically above those calculated from the algebraic fit, shown in Fig.3(b). 
Fits to the exponential form plus its image lead to better results in this case. 
We found that the values of a are quite small and it appears that they become 
smaller as L increases. This suggests that the algebraic contribution will dominate 
in the asymptotic limit. Notice also that 770 increases slowly with L and it does 
not seem to saturate for the larger lattices in both types of fits. The values for 
the exponents calculated from both types of fits in the largest lattices do agree 
with those obtained in Refs. 44 and 46, with the later results given in Table 1C 
for comparison. The value r/o = 0.2386(2) calculated from our algebraic fits is 
in very good agreement with MC renormalization group calculations. 46 However, 
the value 770 = 0.2713(2) obtained by assuming an exponential fit plus its images 



a are smaller than 10 



suggesting a coherence length too large for the fit to be 
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at T = Tbkt also agrees with the results obtained studying the relationship 
between £o(T) and the susceptibility x(T) as T — > T^ KT . The results from the 
present analysis suggest, as expected, that the algebraic form gives a better fit in 
the low temperature phase. 

As a further check to this conclusion, to be used in the FFXYM analysis, we 
now show that the values of Tbkt, Y° and 770, which were determined indepen- 
dently, are consistent with the universal value for the jump in T°(T = Tbkt)- 
We calculated the magnitude of the jump in T°(T = Tbkt), for sizes L = 
8, 16, 24, 32, 48, 60, 72, 84, and 96, and carried out a finite size analysis. We 
used the Tbkt = 0.9035(6) obtained from the high temperature analysis of the 
U(l) correlation functions, to be described below. The results are given in Table 
3A. We got the extrapolated value of 

T°(T B kt) = 0.5986(49) (25) 

for the infinite system, by fitting a straight line to T°(Tbkt) vs L^ 1 , and using 
the L = 48 to 96 results (See Fig. 4(a)). We also estimated the magnitude of the 
jump from the universal intercept of Y°(Tbkt) with the 2Tbkt/k line getting 

T°(Tbkt) = 0.5752(4), (26) 

which is about 3% below the value in Eq(25). Next, we considered the relationship 
between the exponent 770 and the universal prediction for the jump in T°{Tbkt) 38 , 
that is 770 = IvbTbkt /[ZiiT {Tbkt)}- Before inserting the numbers it is important 
to stress the fact that the three quantities appearing in this equation were obtained 
from three different calculations. The critical temperature was calculated from 
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the coherence length analysis in the high temperature phase, T°(Tbkt) from a 
finite size analysis at Tbkt, and the 770 from a finite size analysis of the algebraic 
correlation functions at Tbkt- Plugging in the numbers gives the result 



which indicates that the values of these quantities satisfy within the errors, the 
universality relation for the jump in T° at Tbkt- 

In conclusion, we have shown in this section that our strategy yields reason- 
able quantitative estimates of the critical temperature, critical exponents and the 
magnitude of the jump of T°(Tbkt) in the XYM. It is reassuring that independent 
calculations lead to essentially the same quantitative results. Building from what 
we have learned in this section about the XYM, in the next section we proceed to 
apply the same logic and analysis to the study of the phase transition(s) in the 



4. Critical properties of the fully frustrated 2-D XY model 

In this section we present the bulk of our thermodynamic and correlation func- 
tion results for both the U(l) and Z2 freedoms. We start by discussing the thermo- 
dynamic properties and then we move on to present our results for the correlation 
functions. 




(27) 



FFXYM. 
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4.1. Thermodynamic properties 



(a) U(l) freedoms 

We begin by discussing the helicity modulus T(T). Previous studies 4 ' 39 of 
the FFXYM and the fully frustrated Coulomb gas 17 have indicated the possibility 
that the jump in T(T) may be different from the universal XYM result. To further 
shed light onto this problem we have studied T(T) as a function of temperature for 
different lattice sizes and carried out a finite size analysis of ^u(i) = Y(T = ^U{1))- 
Figure 5 shows the results in the temperature range 0.20 < T < 0.65 obtained from 
runs for L = 8, 16, 32 with 250K MCS and L = 60 with 200K MCS. Notice that at 
low temperatures the finite size effects are almost negligible, however, they become 
important in the critical region. The behavior for L=32 and 60 is about the same 
in the temperature region where T was calculated. To investigate the magnitude 
of the jump in T we proceed as in the XYM calculations. We performed a finite 
size analysis of T at the critical temperature ^\j(i) — 0-44, found from a high 
temperature analysis of the correlations, to be discussed later. The simulations 
were carried out in lattices of size L — 8, 16, 24, 32, 48, 60, 72, 84, and 96 and the 
results are given in Table 3B. The behavior of ^u(i) as a function of 1/L is shown 
in Fig. 4(b) for L = 96, 84, 72, 60 and 48. The value 

T U{1) = 0.37(1) (28) 

was estimated by extrapolating the data to an infinite lattice. This result suggests 
that for the lattice sizes and statistics of our simulations, the jump in the helicity 
modulus for the FFXYM is about 23% below the XYM result. The estimate 
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T(T = Tqg) — 0.34(1) was obtained from MC simulations of the fully frustrated 
Coulomb gas on a square lattice 17 using the formula T = Tgg/I^^cTcg] with the 
values e" 1 = 0.63(3), and Tcg — 0.129(2). Here, e c is the value of the dielectric 
constant at Tcg- Thus, our extrapolated value for T c and the one obtained from 
the Coulomb gas data differ by about 7%, which can be considered in reasonable 
agreement. Note, however, that one cannot rule out the possibility of a smaller 
value of Ttcg f° r larger lattices. Nonetheless, we do not believe that the trend 
would change significantly from the result given here. This result confirms previous 
suggestions 4 ' 17 and gives support to Minnhagen's heuristic conjecture 39 about the 
difference between the jump of Ty CG for the frustrated Coulomb gas in a square 
lattice as compared to the XYM universal jump. 

It has also been suggested 14 that the transition in the FFXYM could be weakly 
first order. To check this possibility we looked at the histogram of the energy about 
and found no evidence for the existence of two competing states. In previous 
MC simulations 4 ' 6 ' 7 ' 14 ' 15 ' 17 it was found that the behavior of the maximum of 
the specific heat as a function of lattice size was consistent with a logarithmic 
divergence, favoring an Ising-like transition. However, MC simulations in larger 
lattices 20 suggest that it is very difficult to distinguish between a logarithmic or 
a power law divergence. We have studied the specific heat and found no signature 
for a logarithmic divergence but we were unable to extract reliable exponents. 

We have also studied the staggered magnetization M s t ag g. (Eq. (5)) due to the 
supercurrents circulating around the plaquettes as a function of temperature and 
lattice size. Figure 6 shows the behavior of M sta gg. as a function of T for L = 16 and 
32. It is non-zero at low temperatures and drops sharply at about T = 0.42. We 
note that finite size effects are almost negligible for these lattice sizes. The behavior 
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of M s tagg. as a function of temperature suggests that it may be considered as an 
order parameter for the U(l) phase transition. Note, however, that the chirality is 
defined in terms of the direction of the circulating currents about the plaquettes 
and thus M s t a gg. can also be thought of as an order parameter for chirality. 

(b) Z2 freedoms 

We calculated the staggered chiral magnetization M. s defined in Eq(7). How- 
ever, A4 S oscillates too irregularly between positive and negative values, thus it was 
more convenient instead to study 7 M 2 S and its fluctuations xt- These quantities 
are plotted as a function of temperature for L=32 and 60 in Fig. 7(a) and 7(b), 
respectively. Ai 2 s goes to unity at low temperatures and it decays sharply to zero 
close to the critical region. Note that x 2 displays an asymmetric behavior close to 
T^ 2 (~ 0.42 for L = 60), where it has a sharp maximum. This indicates that the 
critical exponents for x 2 s above and below Tz 2 should be different. For T > Tz 2 we 
fitted the MC data to 

M S 2 ~(6 Z2 ) 2 ^, and x 2 s ~(ez 2 y 2j (29) 
while for T < Tz 2 X 2 S was fitted to 

X 2 S ~ (-ez 2 )- 2 ^'. (30) 

We extracted the critical exponents 2(3, 27 and 27/ by a straight line fits to \n(M 2 ) 
vs ln(e^ 2 (L)) and In (x 2 ) vs ln(|e^ 2 (L)|) for temperatures within 10% from the 
estimated Tz 2 (L). Here we used the notation ez 2 = (T — Tz 2 )/Tz 2 (L), with Tz 2 (L) 
the temperature at which Ai 2 s goes steeply to zero and x 2 shows a maximum for a 
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given L. The exponents obtained for the largest lattice were 

2/5 = 0.1936(35), 2 7 / = 1.025(79), and 2 7 = 1.82(13) (31). 

These exponents clearly differ from the corresponding 2-D Ising model exponents 
2(3 = 1/4, 2 7 = 2 7 / = 7/2. We note that our chiral order parameter exponent 
does agree with the value 2(3 = 0.20(2) obtained from MC transfer-matrix studies 
of the FFXYM 25 . The results for L = 16, 32 and 60 are given in Table 4. Notice 
the consistency in the behavior of M 2 S which falls off to zero at about the same 
temperature where xi has a maximum, indicating that the Z2 phase transition 
takes place at Tz 2 ~0.42. 

4.2. Correlation functions 

(a) U(l) correlations 

In this subsection we discuss our MC results for the gauge invariant phase cor- 
relation functions obtained from simulations in lattices from L=16 up to L=240, 
with periodic boundary conditions. Some of these results have already been dis- 
cussed in Ref. 22 and thus we will make reference to them here. To reduce finite 
size effects the lattice sizes at each temperature were chosen such that L/£ > 5. 
As we mentioned in Section 3.2 this criterion has proven to work well in the numer- 
ical calculations of correlation functions in the XYM. We showed in Ref. 22 that 
the zero-momentum phase correlation function gu{i){ r ) has an oscillatory behavior 
with period 1/2, which comes from the Aharanov-Bohm phase factors discussed in 
Section 2. 3. 3 At higher temperatures we found that the oscillatory behavior disap- 
pears, as one would expect. As the critical temperature is approached from above 
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the oscillations increase in amplitude and saturate below 7f/(i) ■ This oscillatory 
behavior led us to separate the correlation functions into two components; one 
for the odd and one for the even lattice sites, with their corresponding coherence 
lengths £ and £ e . The MC data for the zero momentum correlation functions was 
fitted to the periodic version of the ansatz given in Eq. (16) for T > This 
procedure incorporates the periodic boundary conditions due to the finiteness of 
the lattice (Eq.(22)). We carried out unconstrained non-linear 3-parameter fits to 
the data to obtain £(T), 77 (T) and the coefficient A for the odd and even correlation 
functions. We followed our XYM approach and fitted the MC data to the linear 
functions 



ln((?(r)) = In A + In 



r 



-ri e -r/Z + (L - r )-'» e -( i -'-)/* 



(32) 



varying i] until a minium for x 2 was reached. In Table I of Ref. 22 we gave the 
results for £ and £ e as a function of temperature and lattice size, as well as the 
statistics of the runs. As one gets closer to the critical temperature from above, 
the coherence length increases exponentially and one needs longer simulations and 
larger lattices in order to get statistically reliable data. Furthermore, the fitting 
parameter rj(T) defined in Eq. (16) increases and oscillates rapidly, for both the 
odd and the even lattices so that an estimate of r\ (Tu^) was not attempted. The 
same situation was encountered in the XYM as described in Section 3.2. We also 
found that as the temperature decreases, A decreases while A e increases, both 
slowly. Far from the critical region we got reliable results for £ and £ e using 
lattices of size L < 60. However, to obtain meaningful results as we got closer 
to the critical region, L was increased keeping the ratio j > 5. For instance, 
we had to increase the size up to L=240 for temperatures that were about 3% 
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and 2% away from T^ufi)- ^ n contrast, in our XYM calculations the approach to 
Tbkt was only on the order of 10%. In Fig. 2 of Ref. 22 we showed the results 
for £ . Similar results are obtained for xi e . We found that for the temperatures 
considered here £ < £ e . In the determination of the critical exponents and the 
critical temperature it is crucial how one fits the data, as discussed in Ref. 46. We 
first tried a 4-parameter unconstrained non-linear fit to the MC data of the BKT 
type (Eq. 17) obtaining the results, 

v e = 0.3133(57), and v = 0.3005(6), (33) 

which are close to 1/3. We then fixed the values v e = v = 1/3, and carried out 
a 3-parameter fit to the data for both lattices. The quality of the fits improved 
and hence we could surmise that the correct value of this exponent may indeed be 
1/3. The first column of Table II in Ref. 22 listed the results obtained from these 
fits together with their corresponding x 2 /dof. For completeness we also carried 
out fits assuming v = 1/2, the standard BKT value, and although the x 2 function 
was smaller than when v — 1/3, we found the differences too small to decide with 
absolute confidence from our data which exponent is the correct one. Nonetheless, 
as it will be seen below, we have other arguments, for example the finite size scaling 
analysis of the data, that yields better results when v — 1/3, suggesting that this 
may very well be the correct value. 

It should be stressed that doing non-linear fits is a non-trivial matter since 
there is no guarantee that the values of the estimated parameters correspond to 
the absolute minimum of the x 2 function. Therefore, one needs to check the results 
very carefully and often times resort to different fitting procedures to cross check 
the results. For instance, a good test to check the stability of the results is to reduce 
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the dimension of the parameter-space by fixing one of the parameters and carrying 
out the fitting procedure for the others. This process should be repeated 43-46 for 
a set of values of the parameter held fixed within an interval about the value which 
is supposed to yield the minimum to x 2 - I n some instances it is worthwhile to use 
the value of more than one parameter, say two, to reduce the problem to a linear 
fit. For example, in the calculation of v and we carried out linear fits to 

HO = HA) + B(T-T u(1) y\ (34) 

with v = 1/2 or 1/3 while varying about the value — 0.446, obtained 

from the nonlinear fits. We found that sometimes these fits led to different values 
of the parameters A and B and also to different values for the minimum of the 
X 2 function. However, the different values of Tj/m extracted from these fits were 
not significantly different. This uncertainty in the analysis must be due to the 
complicated topology of the parameter-space and, although our calculations are 
very extensive, the number of points used in the fits with their corresponding 
statistical significance may not be sufficient to obtain a clear minimum for the x 2 
function. 

Another possible source of problems relates to the size of the errors that weight 
each value of £(T) in the fits. In our calculations most of the errors were on the 
order of 1CT 3 (see Table I of Ref. 22) and, because of the small number (12) of 
available data points we obtained relatively large values for the x 2 function. This 
in turn led to small Q values, (Q being the goodness of fit). In most cases we found 
Q< 0.10, suggesting that the data did not fit the model well. To sort out this 
problem we followed standard practice by setting a\ = 1 and carried out the fits 
again 51 . We found that by doing this the values of Q became larger than 90% in 
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most cases, and the values of the fitted parameters remained basically the same 
but with slightly smaller errors. The values shown in Table II of Ref. 22 were in 
fact obtained following this type of analysis. One will need to have more points 
to improve the quality of the fits. To do this one needs to carry out calculations 
even closer to the critical point. One of the major problems in undertaking such 
a program is the lack of an efficient algorithm that could reduce effectively the 
critical slowing down. For the known algorithms that reduce the critical slowing 
down there has been the need to know in some detail their elementary excitations. 
In a sense this is equivalent to having to know the main features of the solution to 
the model before an appropriate algorithm can be tailored. 

As an additional test of the reliability of the results for Tjj^ and v, we carried 
out a finite size scaling analysis of the data for £ G and £ e . For a finite system, 
assuming periodic boundary conditions, the usual T > Tjj^ finite size scaling 
ansatz for a BKT transition is 

£(T,L) -LF^^exp^e-^j, (35) 

with F(: the scaling function, not known a priori, which must satisfy the conditions 

F^(x) =0, as x — > 0, 

F^(x) < oo, as x — > oo. 

The idea is to find the set of parameters B, v and f° r which the data for 

different temperatures and lattice sizes fall onto one curve. Fixing v = 1/3, we 
varied the values of B and ^[/(i) about their values obtained in the previous fits. 
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We found that as we moved away from those values in the increasing or decreasing 
directions, the data became more scattered. However, very close to the values 
found from the previous fits the data fell very close to a unique curve. The values 
for which the data collapsed approximately onto the universal curve were 

B e = 1.045, B = 0.999, and T£ (1) = 0.440, T£ (1) = 0.442. (36) 

These numbers are in rather good agreement with the values found in the previous 
fits. In the inset of Fig. 2 of Ref. 22 we showed the results of such analyses for the 
odd lattice. Similar results are obtained from the analysis of even lattices. We see 
that close to the critical region the points corresponding to lattices L < 60 are far 
from the universal curve. The equivalent finite size scaling analysis fixing v = 1/2 
always led to a rapidly increasing curve suggesting that closer to the critical point 
it would diverge. This analysis provides further support in favor of v = 1/3. 

As in the XYM analysis, we also tested a power law fit to the £ (T) data, and 
the results are given in the third column of Table II of Ref. 22. We find that the 
BKT and power law fits are of comparable quality, as in the XYM case. Hence, one 
cannot be absolutely sure from this analysis alone which one of the two fits is the 
correct one. In trying to resolve this ambiguity we also calculated gu(i)( r ) below 
mostly for L=60. In fitting the corresponding data to an algebraic form we 
followed a procedure that parallels the one discussed in Section 3.2. The results of 
the analyses are presented in Tables 5 A, 5B. In Figure 8 we show the exponents 
i]o(T) (□) and r) e (T) (o) obtained from the algebraic fits to gu(i)( r )- A careful look 
at the numbers indicates that the trend in r)(T) for both lattices is qualitatively 
similar to the one found in the XYM, but they are quantitatively different. The 
exponential fits to gjj(i)( r ) appear to yield better results with a < 10~ 2 , and 
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larger values for r](T) than the ones obtained with an algebraic fit, see the results 
in Tables 6A and 6B. Nonetheless, the values of a decreased as the lattice sizes 
increased, suggesting that at the asymptotic limit the leading contribution will 
mostly come from the algebraic part of the correlations. We also calculated r) e 
and 7] at the average critical temperature 7}y(i) = 5(^/(1) + ^u(i)) obtained from 
the high temperature analyses for lattices with L = 32, 40, 48, 60, 72, 84 and 96. 
The results from the finite size analysis of the algebraic fits are summarized in 
Tables 5C and 5D, whereas we list the corresponding results for the exponential 
fits in Tables 6C and 6D. The resulting values for r/ e and r\ as a function of L arc 
shown in the inset of Figure 8. For comparison we also show 770. Observe that i] 
is systematically above i] e and that the behavior of these exponents as a function 
of L is qualitatively similar to those found in the XYM, e.g. the r/'s increased 
monotonically with L without appearing to saturate for the values considered. 
However, the ?y's do seem to reach a more asymptotic value for the FFXYM than 
for the XYM. From the above analysis we extracted the results, 

Vo(T u{l) ) = 0.1955(3) and V e(T u(1) ) = 0.1875(3). (37) 

On the other hand we get 

Vo(T {u{1)) ) = 0.2521(3) and V e(T {u{1) ) = 0.2480(3), (38) 

assuming the exponential fits to the correlations. We note that in the algebraic fits 
the values of i] and r\ e are smaller than in the exponential case, as in the XYM 
analysis, and clearly r] ^ 770- 

Again as in Section 2.3, we carried out a check of the universal jump relation- 
ship as applied to the average values of the even and odd lattice results. We found 

33 



that the universal jump relation is indeed satisfied for our FFXYM results with 
the jump at of 

T(T U{1) ) = 0.37(1), (39) 

clearly different than the XYM universal jump. Apart from giving a strong con- 
sistency check of the set of results obtained by independent calculations this is a 
surprising finding, for there is no a priori reason why the universality results should 
be valid in the FFXYM, in particular in light of the non Ising and XY results from 
our study. In the XYM the universality of the jump in T is a consequence of an 
underlying universal RG result 36 ^ 38 . We can then just surmise that there may 
be an underlying RG argument that will lead to an understanding of the physical 
properties of the FFXYM. 

(b) Z2 correlation functions 

Let us now turn to the discussion of the correlation functions for the chiral 
degrees of freedom. Our study here will be less detailed than in the U(l) case, 
mainly concentrating on the temperature region above Tz 2 , although a few results 
for T < Tz 2 will also be discussed. Prior information about the chiral critical 
exponents is available so that we can compare our results to them. The calculation 
of the zero momentum chiral correlation functions defined in Eq. (15) is less 
demanding than in the U(l) case since one expects that £ x diverges algebraically. 
The analysis of g x (r) followed a similar logic to that of the U(l) study. The Fig. 9 
inset shows g x (r) vs r above and below Tz 2 - The results for the coherence length 
£ X (T) for different lattice sizes are given in Table 7. Figure 2 of Ref. 22 showed 
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the results of a power law fit of the data to 



£x ~ (^r^- (4°) 

We also found that inclusion of the errors of g x (r) in the fits yielded results with 
rather low confidence levels (Q < 0.10), as in the analysis of gjj(i)( r )- Again 
we followed standard procedure by assigning the same weight to each data point, 
with the resulting Q values in most cases above 0.90, while the fitted parameters 
remained essentially the same. In the £ x case the errors in the fitting parameters 
became smaller after normalization of the errors in the g x {r) data points. It is 
difficult, however, to be absolutely sure that the values of the parameters found 
correspond to the absolute minimum of the x 2 function, as happened in the U(l) 
case. Therefore, in addition to the nonlinear 3-parameter fits, we also fitted the 
data to the linear function 

\n(^) = \n(A x )-u x \n(T-T Z2 ). (41) 

Using a least square fit and varying the Tz 2 values about 0.42 we found 

A x = 0.33(2) and v x = 0.80(1), for T z , 2 = 0.430, (42) 

with x 2 — 3.99 x 10~ 3 . We also performed least square fits of the data to the 
straight line 

& 1/v *=A x T + b x , (43) 

for fixed values of v x with A x = A x ^ x and b x = A x ^ x Tz 2 - Remarkably, the 
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results obtained were 

A x = 0.36(3), and T Z2 = 0.432(9), for v x = 0.760, (44) 

with x 2 = 9.7 x 10~ 2 . Note that the results in Eq(44) are close to those in Eq(42), 
with essentially the same Tz 2 - In Table 8 we give the values of the parameters 
extracted from the nonlinear fit. Our result for v x agrees quite well with recent 
finite size scaling analysis 20 that gave v x = 0.85(3), as well as with the MC transfer 
matrix calculations 25 . The advantage of the finite size scaling analysis is that v x 
was obtained from a one parameter fit without needing a precise value for T^ 2 , as 
in our analysis. Therefore it appears that the v x and Tz 2 values obtained here from 
the nonlinear fits may in fact be very close to the correct ones. It is important 
to emphasize that the Tz 2 found here is consistent with the temperature at which 
M? s fell to zero, and Xs displayed a sharp maximum. 

In summary, our numerical analysis of the chiral degrees of freedom led to the 
critical exponents 

2/5 = 0.1936(35), 27/ = 1.025(79), 2 7 = 1.82(13), and v x = 0.875(35). (45) 

These results strongly indicate that the Z2 phase transition is not an Ising-like 
transition as had been suspected from previous thermodynamics studies of this 
model. Note that in our calculations the difference between the Tz 2 and Tjjn\ : is 
about 7%, which may not be considered as different within the size of our estimated 
errors. Equivalently, one cannot rule out the possibility that in improved numerical 
simulations and closer to the critical point this difference may disappear. 
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5. Conclusions and outlook 



In this paper we have presented results from extensive MC calculations of the 
FFXYM. We have explicitly analyzed the separate contributions from the U(l) and 
Z2 freedoms. We have extracted the U(l) and Z2 critical exponents from direct 
calculations of their corresponding correlation functions and selected thermody- 
namic properties. We found compelling quantitative evidence that the U(l) and 
Z2 critical exponents are clearly different from those of the usual 2-D XY and Ising 
models. We have tested our results using several consistency checks. Our result for 
the Z2 correlation length exponent v x is essentially the same as the one obtained 
from other independent numerical calculations 20 ' 25 . There are no previous calcu- 
lations of the U(l) exponents to with which to compare our results. However, a 
reanalysis of the experimental data leads to an 77 exponent that is clearly different 
from the XYM result and that agrees reasonably well with the one found in our 
calculations. 

Our results strongly suggest nontrivial critical behavior in the FFXYM, in 
which the U(l) and Z2 freedoms are coupled in a way so as to yield novel critical 
exponents. We leave for the future the question of producing the physical under- 
standing of the results presented in this paper, in particular the apparent relation 
between the U(l) results, ?7o _1 (/ = 0) = 4, ^o _1 (/ = 0) = 2 and our values, 
rj-^f = 1/2) = m ~ l {f = 0) + 1 and v~\f = 1/2) = VQ~\f = 0) + 1, together 
with the validity of a universal jump for the / = 1/2 helicity modulus. 

In spite of the extensive calculations and detailed analyses carried out in this 
paper, improved MC simulations of this system need to be done. More data at 
temperatures closer to the critical point are required to be able to obtain more 
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accurate estimates of the critical exponents and critical temperatures. As noted 
above, among the major limitation in this program is the lack of a MC algorithm 
that could effectively minimize the critical slowing down. Another important limi- 
tation is that calculations of correlation functions near criticality require ever larger 
lattices which sharply increases the computer power requirements. 
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Figure captions 



Figure 1. Possible trajectories r^(r, r t) and r#(r, f /) used in evaluating 
9u(i)(r,rt,r A ) and ^(i)(f, f /, F B ). 

Figure 2. Correlation length £o as a function of temperature. The circles 
denote the results from the fits to the correlation function given in Eq.(16), with 
their corresponding statistical errors. The continuous line is the result of a fit of 
the data to the expression for £o(T) given in Eq(17). The specific values of the 
fitting parameters is given in Table 1A. The dashed line was obtained using the 
fitted parameters from Ref. 43, 

Figure 3. (a) Results for rjo(T) obtained from algebraic fits to go(r) for L=60. 
The dotted line is the spin-wave results, (b) shows i]q c (T = Tbkt) versus L. 

Figure 4. (a) Finite size analysis of T°(T = Tbkt) = ^°t B kt-> w ^ Tbkt 
obtained from the go(r) analysis, as a function of 1/L for L = 96, 84, 72, 60, 48, 
and 32. The straight line is a linear fit to the data, with the L = oo extrapolated 
value indicated, (b) The same as in (a) for T(T = Tjj^) = Tt u{1) - 

Figure 5. T as a function of T for different L sizes. Note that the data for 
L = 32 and 60 almost fall on top of each other, suggesting that the L dependence 
is almost negligible for L > 32. 

Figure 6. Staggered magnetization due to the superconducting currents defined 
in Eq. (5) as a function of T for L = 32 and 60. The fall off to zero occurs at 
about T = 0.42. 

43 



Figure 7. Staggered chiral magnetization square M 2 (a) and susceptibility Xs 
(b) as a function of T for L = 32 and 60. The results for both lattice sizes are 
essentially on top of each other. Note that M 2 S goes to zero at essentially the same 
T at which xi has a maximum with Tz 2 ~ 0.42. Observe the asymmetric behavior 
of Xs about critical temperature. 

Figure 8. Results for r] (T) (□) and rj e (T) (o) obtained from algebraic fits to 
9u(i)( r ) f° r L=60. The inset shows the results for the finite size analysis for the 
critical exponents t]q (□), if (o) and i]° (x), at criticality. 

Figure 9. Data for the chiral coherence length £ X (T) (x) calculated from g x (r). 
The solid line is the fit to the form given in Eq(20). The inset shows g x (r) as a 
function of r for temperatures above and below Tz 2 . 
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Table captions 



Table 1. Critical exponents and critical temperatures for the XYM. (A) In the 
first and second columns we give the results obtained from fitting the T > Tbkt 
£o data to Eq(17) and Eq(20), respectively. The third column gives the parameters 
obtained from an algebraic fit to go(r) at T = Tbkt- (B) Same as in (A) for the first 
two columns with results from Ref. 46, obtained from unconstrained 4-parameter 
nonlinear fits. (C) The first two lines are the t]q exponents obtained from a high 
temperature analysis of the susceptibility x (from Ref. 44). The third and fourth 
lines give t]q obtained from MC renormalization group calculations (Ref. 46). 

Table 2. (A) Results from fits to the data for the correlation function go(r) to 
the form given in Eq.(18) for L=60. (B) L dependence of 770 and Co at T = Tbkt- 
(C) Results from exponential fits to the data of go(r) for T < Tbkt f° r L=60. (D) 
L dependent results for i]q, ao and Co from exponential fits to the correlations at 
T = T B kt- 

Table 3. (A) Gives the finite size results for Y° BKT , while (B) gives the 
corresponding results for Ty . 

Table 4. Magnetic critical exponents for the Z2 transition for L = 16, 32 and 
60. See text for definition of the parameters. 

Table 5. Results from algebraic fits to <7j/(i)(r), (A) odd and (B) even lattices 
for T < Tf/fi) with L = 60. Results of a finite size analysis of the algebraic fits to 
9u(i)( r ) f° r ^e odd (C) and even (D) lattices at T = Tjj^y 

Table 6. Results from exponential fits to gu(i){ r ) f° r the odd (A) and even 
(B) lattices for T < with L = 60. Results from a finite size analysis of 
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exponential fits to gu(i)( r ) f° r the odd (C) and even (D) lattices at T = 

Table 7. Results for £ X (T) obtained from g x (r) at different temperatures and 
lattice sizes. 

Table 8. Exponent v x and critical temperature Tz 2 obtained from nonlinear 
fits to a power law divergence. 
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